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Slowing down atomic diffusion in subdwarf B stars: 
mass loss or turbulence? 



o 

(N 
< 

in 



6 

in 
c3 



> 
00 

m 

06 
o 



Haili Hu 1 *, C. A. Tout 1 , E. Glebbeek 2 and M.-A. Dupret 3 

1 Institute of Astronomy, The Observatories, Madingley Road, Cambridge CB3 OH A 

2 Department of Physics & Astronomy, McMaster University, 1280 Main Street West, Hamilton, Ontario, Canada L8S J { M1 
3 Institute d'Astrophysique et Geophysique de I'Universite de Liege, Allee du 6 Aout 17, 4000 Liege, Belgium 



Accepted 2011 July 20. Received 2011 July 19; in original form 2011 June 14 



ABSTRACT 

Subdwarf B stars show chemical peculiarities that cannot be explained by diffusion 
theory alone. Both mass loss and turbulence have been invoked to slow down atomic 
diffusion in order to match observed abundances. The fact that some sdB stars show 
pulsations gives upper limits on the amount of mass loss and turbulent mixing allowed. 
Consequently, non-adiabatic asteroseismology has the potential to decide which pro- 
cess is responsible for the abundance anomalies. We compute for the first time seismic 
properties of sdB models with atomic diffusion included consistently during the stellar 
evolution. The diffusion equations with radiative forces are solved for H, He, C, N, 
O, Ne, Mg, Fe and Ni. We examine the effects of various mass-loss rates and mixed 
surface masses on the abundances and mode stability. It is shown that the mass-loss 
rates needed to simulate the observed He abundances (10 -14 < M/M yr _1 < 1CP 13 ) 
are not consistent with observed pulsations. We find that for pulsations to be driven 
the rates should be M < 10 -15 M Q yr _1 . On the other hand, weak turbulent mixing of 
the outer 1CT 6 M can explain the He abundance anomalies while still allowing pul- 
sations to be driven. The origin of the turbulence remains unknown but the presence 
of pulsations gives tight constraints on the underlying turbulence model. 
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1 INTRODUCTION 



X 



Subdwarf B (sdB) stars are subluminous B-type stars that 
show strong broad Balmer lines and weak or no He I absorp- 
tion lines. They are ubiquitous not only in our own Galaxy, 
where they are found in all stellar populations, but also in 
old galaxies where they are believed to cause the UV-upturn 
( Brown et al.|1997| ). [Heber| ( |1986| linked the sdB evolution- 
ary stage to that of Extreme Horizontal Branch (EHB) stars 
which are core He burning stars with an unusual thin H en- 
velope (A/ cnv < 0.02 Mq). Such a star must have formed 
either from a red giant that lost nearly its entire H envelope 
or from a He white dwarf merger. Although the formation 



2003 and references 



channels are well studied (Han et al 
therein), they cannot explain all observed features such as 
binary period distribution and slow rotation. A promising 
way to discriminate between different formation channels is 
with asteroseismology ( Hu et aL|[2~0 08) but, before this can 
be achieved, the theoretical models need improving. Cur- 
rently, the uncertainties in the models are larger than the 
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observational errors obtained with space missions such as 
CoRoT and Kepler. This paper is part of our ongoing work 



(Hu et al. 2009 20101 to improve the input physics used in 



seismic modelling of sdB stars. 

Interestingly, two subgroups of sdB stars are pulsating 
with variable class names V361 Hya and V1093 Her. The dis- 
covery of short-period (100 — 400 s) pulsations in V361 Hya 



stars by Kilkenny et al. 



( 1997 \ coincided with the predic- 



tion of unstable pressure(p)-mode pulsations in sdB stars by 



Charpinet et al. ( 1996 1. The excitation was attributed to the 



opacity mechanism enabled by Fe accumulating diffusively 
in the stellar envel ope (|Charpinet et aL||1997| ). A few years 
later, Green et al. ( 2003 1 discovered long-period pulsations 
(30 — 120 min) in the cooler V1093 Her stars. Fontaine et al. 
( |2003[ ) showed that the same opacity mechanism could also 
excite long-period gravity(<7)-modes. However, they found 
too cool a theoretical blue-edge of the (/-mode's instabil- 
ity strip compared with observations. The problem can be 



partly solved by using OP opacities ( Badnell et al. 2005 1 



instead of OPAL's (Iglesias & Rogers 



ing that Ni accumulates as well as Fe ( Jeffery & Saio 2006 



1996) and by assum- 



Furthermore, inclusion of H-He diffusion shifts the theoret- 
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ical blue-edge even closer to the observed value of 30 kK 
( Hu et al~||2009 1. It is evident that correct modelling of the 
V1093Her instability strip awaits seismic models with time- 
dependent diffusion such as we produce here. 

Also interesting is that sdB stars are chemically peculiar 
without significant spectroscopic differences between vari- 
able and constant stars ( jO'Toole fc Heber|20 06 Blanchcttc 
et al. 2008). Typically He is deficient with number ratios 



10 nHe/«-H ^ 0.1 and there is a trend for He abun- 

dance to increase with T eff ( |Edelmann et al.||20"03"| ). The Fe 
abundance is close to solar irrespective of stellar popula- 
tion and T ff whereas other iron-group elements, such as Ni, 
can be strongly enriched. Light metals (C, N, O, Ne, Mg) 
are usually depleted, although there is some scatter between 
stars (Edclman n et al.pO OG Gcicr ct al. 2010 1. The observed 
abundance anomalies can be explained in part by atomic 
diffusion acting in the radiative envelope. Atomic diffusion 
alone, however, would cause all surface He to sink on a very 
short time-scale (about 10 4 yr) compared to the EHB evo- 
lutionary time-scale (about 10 8 yr). Hence there must be a 
competing process which is generally accepted to be mass 
loss. Recently turbulence has been proposed to play this role 
instead ( |Michaud et al.|2011 l. 

The observed He abundances have been reproduced 
with models that include atomi c diffusion and mass-loss 
rates of 10" 14 ln " 13 



10"" M yr" 1 (Fontaine fc Chayer 



1997 



Unglaub fc Bues|2001 1. For higher rates the effects of diffu- 



sion become unnoticeable, whereas for lower rates He would 



sink too quickly. Vink & Cassisi ( 2002 1 presented a mass- 



loss recipe for (E)HB stars derived from line-driven wind 
models. Their results give upper limits for the mass loss 
rates (M < 10" 11 Mg yr" 1 ) and might apply to the most 
luminous sdB stars that show anomalous Ha line profiles 
( Heber et aLp003"l |Vmk|[2004l ) . However, for most sdBs no 
observational evidence of mass loss has been found so far. An 



independent wind model by Unglaub ( 2008 1 showed that for 
weak winds, M < 10" 12 Mg yr" , the metals decouple from 
H and He and, for rates below 10" 16 Mg yr" 1 , the winds be- 
come purely metallic. Such selective winds could be respon- 
sible for some of the observed metal abundance anomalies 
but they cannot explain the observed He abundances. 



Alternatively, Michaud et al. (20111 showed that tur- 



bulent mixing of the outer 10" Mg could also produce 
most of the observed abundance anomalies. Unfortunately, 
their computational setup restricted their calculations to 
the first 32 Myr of the E(HB) and to He abundances above 
A (He) w 10~ 4 . They do not give a physical origin for the 
turbulence but simply mix the outer layer in order to repro- 
duce the solar Fe abundance observed in most sdB stars. In 
their models Fe reaches solar abundance in the entire mixed 
region. This could prevent the driving of pulsation modes. 
Also mass loss works against the driving mechanism. 



According to Chayer et al. (20041) and |Fontaine et al. (20061 
mass- loss rates above M w 6 x 10" 15 Mg yr" would grad- 
ually destroy the Fe reservoir and an sdB pulsator becomes 
constant within 10 7 yr as it evolves. Their results are based 
on a non-diffusive stellar model assuming an initial Fe pro- 
file from an equilibrium between gravitational settling and 
radiative levitation. The question is, however, can Fe build 
up in the first place in the presence of stellar winds? If so, 
shouldn't diffusion take place at the same time that mass 
is lost and help build up the Fe reservoir again? To answer 



these questions we perform a consistent analysis where we 
compute the effects of mass loss and atomic diffusion simul- 
taneously. Furthermore, we determine which process (mass 
loss or turbulence) is responsible for retarding atomic diffu- 
sion in sdB stars by performing a pulsational stability anal- 
ysis. Since non-adiabatic effects are responsible for mode 
driving and damping, such a study is termed non-adiabatic 
asteroseismology. 

Beware that the term 'atomic diffusion' is not always 
used consistently in the literature. For clarification, in this 
work we include all physical processes that contribute to 
atomic diffusion. These are gravitational settling, thermal 
diffusion, concentration diffusion and radiative levitation. 
The corresponding abundance changes are caused by pres- 
sure gradients, temperature gradients, concentration gra- 
dients and radiative forces, respectively. Indeed, we treat 
atomic diffusion in a multi-component fluid ( Burgers|1969 1, 
using diffusion coefficients derived from a screened Coulomb 
potential ( Paquette et al.|1986 l. We account for partial ion- 
ization by using a mean ion charge per element. We calculate 
radiative forces from first principles using atomic data. Our 
approach gives a high accuracy treatment of atomic diffu- 
sion because it is free of many common assumptions, such 
as full ionization, simplified treatment of trace elements and 
approximate radiative forces. 

We explain our method in Sect. [2] In particular the 
computations of stellar evolution (Sect. |2.1[ ), radiative ac- 
celerations (Sect. 2.2 1 and atomic diffusion (Sect. 2.31 are 
described. In Sect. |3| we present the results for typical sdB 
models. The effects of atomic diffusion (Sect. 3.1 1, mass loss 
(Sect. 3.2 1 and turbulence (Sect. 3.3 1 on the abundances and 



mode stability are evaluated. We summarize and discuss our 
main results in Sect. [4] We also give a brief discussion of the 
possible physical origin for turbulence. 



2 COMPUTATIONAL METHOD 

This work presents the first seismic computations of stellar 
models with a self-consistent treatment of atomic diffusion. 
The changes in the chemical abundances are accounted for 
by continuously updating the radiative accelerations and the 
Rosseland mean opacity, not only during the stellar evolu- 
tion but also in the pulsation calculations. We describe here 
the computational tools and methods we use to accomplish 
this. 



2.1 The stellar evolution code 

We compute stellar models with the stellar evolution code 
STARS ( |Eggleton|1971[). T his code has been frequently up- 
dated (e.g. |Pols et al.|1995| ) and many derivatives are circu- 
lating. In our version the necessary modifications have been 
made to construct stellar models suitable for seismic stud- 
ies with the non-adiabatic pulsation code MAD written by 
Dupretj pOOl]) (see |Hu et "aJT||2008[>. Furthermore, the im- 



plementation of gravitational settling, thermal diffusion and 
concentration diffusion is described in |Hu et al.| ( [2010| . In 
the present work we add the process of radiative levitation. 
We follow the abundance changes of H, He, C, N, O, Ne, 
Mg, Fe and Ni while the remaining minor species are kept 
constant. 
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2.1.1 Input physics 



We take nuclear reaction rates from Angulo & et al. ( 1999 1, 
except for the 14 N(p, 7) ls O rate for which we use the rec- 
ommended value by |Herwig et aT] ( |2006[ ) and |Formicoia| 
fc LUNA Collaboration! " 2002 1. Neutrino loss rates are ac- 
cording to Itoh et al. (|1989 19921. Convection is treated 



with a standard mixing-length prescription ( |Bohm-Vitense 



1958) with a mixing length to pressure scale height ratio 



of I /Hp — 2.0. Convective mixing is treated as a diffusive 
process in the framework of mixing length theory. 

The Rosseland mean opacity kr, is computed in-line 
during the evolution to account for composition changes con- 
sistently. For this we use data and codes from the Opacity 
Project (OP, |Seaton et aL]|1994| and |Badnell et~aT1|2005 l 
which also allows computation of radiative accelerations g la d 



(see Section 2.2|. This is in contrast to all previous work 
with the STARS code that made use of interpolation in pre- 
built opacity tables. For high temperatures outsi de the OP 
range (T > 10 8 K), we still use OPAL tables (Iglesias fc 



Rogers 1996 I combined with conductive opacities ( |Cassisi 
et al. 12007) 



2.1.2 Simultaneous solution 

The STARS code distinguishes itself from other evolution 
codes by its fully simultaneous, implicit and adaptive ap- 
proach. In other words, the equations of stellar structure, 
composition and mesh are solved together, at each timestep 
and each iteration, and variables are from the current 
timestep. Within the scheme of an adaptive mesh, mass loss 
is simply included by adapting an outer boundary condition. 
More importantly, we find that such an approach is free from 
some of the notorious numerical instabilities that accom- 
pany diffusion calculations. Such instabilities occur because 
atomic diffusion, and in particular radiative levitation, can 
happen on a much shorter time-scale than the evolution of 
the stellar structure, at least in the outermost layers. In non- 
simultaneous approaches, where the composition is solved 
separately from the structure, this can be dealt with by al- 
lowing multiple composition timesteps within one evolution 
timestep (see e.g. |Turcotte et al.|[l998[ ). In addition, evo- 
lution timesteps should remain short if diffusion velocities 
are kept constant during a timestep. This requires lengthy 
calculations because the radiative accelerations are compu- 
tationally expensive to evaluate. 

In our calculations, the timesteps are short (10 2 yr) at 
the start of the evolution but can grow to 10 yr, which is of 
the same order as for models without diffusion. The timestep 
size is determined by the change in variables from one 
timestep to the next. Normally, the temperature and degen- 
eracy have the largest influence but if diffusion acts rapidly 
the abundance variations limit the timestep size. One might 
wonder whether the results are sensitive to the timestep size, 
particularly if mass loss is included. We checked, for one 
simulation with mass loss, that the results do not change 
noticeably if timesteps are kept below 10 3 yr. 

Unfortunately, our numerical scheme requires some ef- 
fort as well. Like most variables, the radiative accelerations 
are calculated implicitly; they are not kept constant during 
a timestep but can change from one iteration to the next. 
In order to find the correct stellar model, we need to know 



how a change in p ra< j affects the stellar structure. In other 
words, we need the derivatives of each g r ad with respect to 
the variables that represent the stellar model. STARS nor- 
mally computes the derivatives numerically but this gives 
instabilities in the case of g ra d, probably because of their ir- 
regularity. Hence their derivatives (or at least some of them) 
must be obtained analytically. This is done by the OP rou- 
tines with some modifications, see Sect. 12.21 



2.2 Radiative accelerations 

We compute radiative accelerations g Ta d and Rosseland 
mean opacities kr using OP atomic data. The OP project 



also provides a set of codes called OPserver (Seaton 2005 



Mendoza et al. 20071 which we have rewritten to suit our 



specific needs. The most relevant modifications are summa- 
rized here. 



2.2.1 Radiative accelerations in stellar interiors 

For Rosseland mean optical depths tr ^ 1, the radiative 
acceleration for element fc can be expressed as follows 

M K R ' /-t \ 

5rad,fe = -. ~7*> I 1 ) 

HkC A-Kr 2 - 

where [i is the mean atomic weight, fit is the atomic weight 
of element fc, c is the speed of light, kr is the Rosseland 
mean opacity, I is the luminosity and r is the radius. The 
dimensionless parameter yk depends on the monochromatic 
opacity data by 

(o>(m)[1 - exp(-w)] - a k (u))du 



7fc 



where a k is the cross-section for absorption or scattering of 
radiation by element fc, a(k) corrects for electron scattering 
and momentum transfer to electrons and fk is the num- 
ber fraction. Eq. |l| is equivalent to equation (6) of Seaton 
( 2005 ) but makes use of the luminosity rather than the effec- 
tive temperature and stellar radius. The latter two surface 
quantities are not known during the iterations because the 
evolution code solves the stellar structure from the centre to 
the surface. 



2.2.2 Radiative accelerations for multiple elements 

For each call OPserver allows the computation of gr ra d for one 
selected element in the mixture and for multiple mixtures 
with varying abundances of fc. Because we are interested in 
the diffusion of nine elements it saves CPU time to compute 
<?radS for multiple elements but only for one mixture in each 
call. 



2.2.3 Mean ion charges 

As mentioned before, elements are treated as if they have an 



mean ion charge Z. It has been said by 



Michaud & Richer 



( 2008 ) that a disadvantage of OP is that the database does 
not contain mean charges which are needed to determine the 
diffusion velocities. We found, however, that Z can easily be 
calculated from their data. That should not be too surprising 
because the level populations of the ionization states are also 
needed to compute g ra d- 
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2.2.4 Interpolation method 

To obtain g ra d and kr for a specified temperature T and 
density p (converted to electron density 7V e ) OPserver inter- 
polates within a tabulated mesh of (T, 7V e ) values. OPserver 
uses two different interpolation schemes. The first is a refined 
bi-cubic spline interpolation ( Seaton|1993 \ within the entire 
(T, N e ) mesh. This scheme is applicable for a single chemical 
mixture only. Secondly, when the chemical abundances vary 
with stellar depth, a less refined bi-cubic interpolation oc- 
curs within a two-dimensional array of (4, 4) points around 
the specified (T, p) values. It should be clear that the second 
scheme applies if we want to take into account the variation 
of abundances owing to atomic diffusion. Unfortunately, this 
scheme is too crude for pulsation studies for which it is re- 
quired to have smooth opacity derivatives 



d log KR 
aiogT 



and 



d log KR 
dlogp ' 



Therefore, we use the general idea of the second scheme 
but replace the interpolation method by that developed by 
Dupret ( |2003[ l for the purpose of non-adiabatic pulsation 
computations. This is a local interpolation method by a 
polynomial of degree 3. It ensures continuity of kt and n p 
and their derivatives. Note that, in our scheme, interpola- 
tion is made within an array of dimension (6, 6) rather than 
(4,4). This increases CPU time, so we only use the third 
interpolation scheme for the pulsation calculations. For the 
evolution calculations OPserver's second scheme suffices. To 
demonstrate the differences between the three interpolation 
methods, we perform test runs on the same stellar model as 
used by Seaton| ( |2005 1 . While kr (not shown) itself is not 
visibly affected by the interpolation method, its derivatives 
(shown in Fig. [TJ are. 



2. 2. 5 Computation of derivatives 

As explained in Sect. |2.l] analytical derivatives of each <? ra d 
are needed in order for the evolution code to converge to- 
wards the correct stellar model. It improves the numerical 
stability to do this for kr and Zi as well. The derivatives we 
compute are 



8 log ffrad.z d log ff r . 

dlogT 



d log ffr; 



aiogp dXi 

d log kr d log kr Okr 
aiogT ' aiogp ' dXi' 
dZi , dZ, 



for 



dlogT 



and 



d log 



The opacity derivatives are already output by OPserver 



(see Sect. 2.2.4 1. The others we added. Strictly speaking, the 
derivatives with respect to logp and logT are not analytic. 
They are however obtained from the interpolaton curve and 
so have an analytic form. 




6.4 6.6 6.8 

log (T/K) 

Figure 1. The opacity derivatives k p {left y-axis) and {right 
y-axis) throughout the star as a function of temperature. OP1 and 
OP2 indicate OPserver's first and second interpolation scheme, 
respectively. 



2.3 The diffusion equations 

Having obtained the radiative accelerations <? ra d, we put 
them in Burgers' diffusion equations (Burgers 1969), 

dpi 



dr 



+ Pi (d — ffrad.i) — UiZieE 



V Kij {wj - Wi) + V K ijZij ra^' 

including the heat flow equations, 
. N 



(2) 



{mi + rrij) 2 

N 



*j ; V ^3 

m,i + mj 



- S~] y— "'^ ^ 2 {imi + rrijZij + Q.8mimjz"j)ri 



+E (mi+mj)2 ( 3 +^~ - 8 ^ 



(3) 



In addition, we have two constraints, current neutrality, 

^ ZiUiWi = (4) 

i 

and local mass conservation, 

^^miUiWi = 0. (5) 

i 

In the above 2N + 2 equations the quantities pi, pi, rii, 
Zi and mi are the partial pressure, mass density, number 
density, mean charge and mass for species i, respectively. 
The total number of species (including electrons) is TV. The 
IN + 2 unknown variables are the N diffusion velocities Wi, 
the N heat fluxes rt , the gravitational acceleration g and the 
electric field E. The resistance coefficients Kij, Zij, z[j and 



Zij are taken from Paquette et al. ( 1986 1. We solve Eqs \2 



( 5[) wi th an adapted version of the routine by Thoul et al 
( 1994] ). In Appendix [A] we describe how we modify Thoul 



et al.'s routine to include radiative forces and Paquette et 
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al.'s resistance coefficients. The output of the routine are 
the diffusion velocities Wi. These are then inserted as an 
advection term in the diffusion equation governing the time 
evolution of abundance Xi , 



dXi 



dt 



+ Ri 



J_d_ 

r 2 Q r 



xdXj 
dr 



dr 



(pr XiWi 



(D + D T )pr- 
pr^ Or v Or J pr^ 

where Ri is the change owing to nuclear reactions, D is the 
combined diffusion coefficient for convection, overshooting 
and semi-convection ( |Eggleton|1972 1 and Dt is the diffusion 
coefficient for turbulent mixing. We take Dt from |Michaud| 
et al.l d2011b, 



where 



J D(He) = [ 



3.3 x 10" 15 T 2 



(6) 



4p In (1 + 1.125 x 10" 16 T 3 p) 



is an approximation to the He diffusion coefficient at a cer- 
tain reference depth. The subscript indicates the reference 
depth in terms of the outer mass coordinate AM = M*—M r . 
The assumed Dt is C ti mes the He diffusion c oefficient at 
AMo and varying as p~™. Michaud et al. (2011 1 used a tur- 
bulence model with C = 10* , n = 4 and AM = 10" 7 ' 5 M Q 
but we shall vary these parameters to examine the influ- 
ence of the efficiency of turbulent mixing on the model (see 
Sect. HO I. 
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Figure 2. Radiative accelerations in the stellar envelope as 
a function of mass coordinate (bottom axis) and temperature 
(top axis). Gravitational acceleration is represented by the solid 
line. This model is at the ZAEHB with T cff = 30, 000 K and 
log(g/cm s~ 2 ) = 5.8 at the surface. 



3 RESULTS 

We investigate the effects of atomic diffusion, mass loss, 
and turbulence in a typical sdB star with total mass M, = 
0.46 M and envelope mass M onv = 3.5 x 10~ 4 M . M cnv is 
the amount of mass above the core boundary and is defined 
at ZAEHB only, since the H profile can change with time. 
At the core boundary X(H) starts to increase from to 0.7 
at the surface. The initial H profile is fitted to the profile on 
the RGB according to |Hu et al.1p009| . The ZAEHB model 
has a T c ff = 30 kK and log(g/cms _ ) = 5.8, and it has a 
metallicity at all depths of Z = 0.02 with metal mixture of 
Grevesse & Noels ( 1993[ ). Although sdB stars could have a 
range of original metallicities, Michaud et al. (2011 i's dif- 



fusion calculations showed that the surface abundances are 
poorly sensitive to the initial metallicity. Thus we do not ex- 
pect our results to be much affected by the choice of initial 
abundances. 

Starting at the ZAEHB we compute the stellar evolu- 
tion until the end of core He burning under different as- 
sumptions for mass loss and turbulent mixing. Thereafter, 
we compute non-adiabatic oscillations of the sdB models. 
We restrict ourselves to modes of spherical degree I ^ 2, 
because higher I values are geometrically disfavoured for 
observation. We search for eigenfrequencies in the range of 
0.1 ^ tj ^ 20 where u is the dimensionless angular eigen- 
frequency, uj — 27r/rd yn , / is the frequency and Td yn is the 
dynamical time-scale. 

3.1 Atomic diffusion, no mass loss/turbulence 

In Fig. [2] we show radiative accelerations for C, N, O, Ne, 
Mg, Fe and Ni in the ZAEHB model. H and He are not 



shown because their <? ra d are so low that they lie outside the 
range of the plot. Because gravitational settling and radia- 
tive levitation are the dominant diffusion processes, gravity 
is also plotted for comparison. 

Fig. [2] tells us that Fe and Ni are levitated, whereas 
the other elements sink sooner or later. This is confirmed 
in Fig. [3] which illustrates the time evolution of the interior 
abundances from ZAEHB to 10 8 yr. Notice that Ni accu- 
mulates to mass fractions comparable to (or even higher 
than) Fe in the outer layer, while it starts out with a much 
smaller mass fraction. This implies that Ni is (at least) as 
important for mode driving as Fe. Notice also the develop- 
ment of an iron-group convection zone within 10 4 yr around 
AM = 10 -10 Mq. The abundance profiles in this region are 
flattened and do not have a detailed shape as predicted by 
equilibrium diffusion theory. 

In the upper panel of Fig. |4j we show the computed 
pulsation periods as a function of stellar age (bottom axis) 
and T c g (top axis). We also plot the ranges of observed peri- 
odicities. We see that this particular simulation represents a 
hybrid pulsator, showing both unstable p- and g-modes. The 
total period range of unstable modes matches the observed 
range (100 s - 2hr) well. In the transition region between 
the p- and gr-modes (400 s - 30min) the agreement is not 
as good, but unconfirmed low-amplitude peaks have been 



observed in that region (Baran et al. 2011 and references 
therein). 

It is important to realize that the low-degree (£ = 1,2) 
g-modes would not be driven at such high T c g ~ 29 kK 
without Ni diffusion. We find that Ni plays an even bigger 
role in mode driving than first suggested by Jcffery & Saio 
(2006). In their models only £ > 2 modes, which are unlikely 
to be observed, are driven at such high T^gs. This is because 
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Figure 3. Interior abundance profiles at different stellar ages. The horizontal dashed lines give initial abundances at ZAEHB, and the 
solid curves are at EHB ages log(f/yr) = 4, 5, 6, 7 and 8 as indicated in the legend. Note that in the case of hydrogen we plot 101ogX(H) 
for visibility. These simulations are with neither mass loss nor turbulence. 



Ni accumulates more than Fe in the driving region which was 
not realized before without diffusion calculations. We shall 
present a detailed investigation of the sdB instability strips 
elsewhere (Hu et al. in prep.). 



The lower panel of Fig. [4] shows the time evolution of 
the surface abundances. The problem of the He abundance is 
clearly illustrated because riHe/nH decreases to below 1CP 4 
within 10 4 yr, or 0.01 per cent of the EHB lifetime, w hereas 

2 (see e.g 



Edel- 



the average observed value is around 10 
mann et al.||2003 l. The other surface abundances are quali- 
tatively in agreement with observations. C, N, O and Fe are 



within the observed ranges given by Geier et al. (20101. Ne 



and Mg are almost completely depleted in our models but 
this could be consistent with measurements that give only 
upper limits. Ni is 1.5 dex above solar. This is on the high 
side but still in agreement with observations that go up to 
1 - 2 dex above solar ( |Q'Toole fc HeberpiOo) . For the re- 
mainder of this work, we shall discuss only the abundances 
of He, Fe and Ni, since only these are relevant for our ar- 
gumentation for constraining the mass loss and turbulent 
mixing. 



3.2 Atomic diffusion and mass loss 

The fact that we observe sdB pulsations puts an upper limit 
on the mass-loss rates, at least for the pulsators. After all, 
Fig. [4] shows that it takes over 10 4 yr before enough iron- 
group elements have accumulated to drive pulsations. The 
accumulation at this age goes down to AM ~ 10~ 10 Mq in 
our models (see Fig. [3JI. So if this outer mass is removed 
within 10 4 yr, or equivalently M > 10 -14 Mq yr -1 , then 
there is no chance for the metals to build up. 

To verify this estimate we perform simulations with 
M = 10~ 15 , 5 x 10~ 15 , 10~ 14 and 10~ 13 M yr' 1 . We as- 
sume that the mass loss is constant and chemical homo- 
geneous. The possibilities of a variable and selective mass 
loss are discussed later. We found that numerical instabili- 
ties occur for models that develop an iron-group convection 
zone and experienc e mass loss at the same time (see also 
Charbonneaull 1993k. This happens for M = 10~ 15 M yr" 1 



in our models. To reduce instabilities we assume, only in 
this case, overshooting of the iron-group convection zone to 
the atmosphere. The extra mixing occurs only in the outer 
10" 11 Mq and has a negligible effect. Still, this simulation 



Slowing down atomic diffusion in subdwarf B stars 7 



29.4 
4 



29.3 



T eff /kK 
29.1 28.9 



28.5 



27.8 



Unglaub & Bues 2001 1. This can also be seen in Fig. [5] 



3.5 



2.5 



5 -3 



II 



iiiiiiiiiiiiiiiiiiii 5 illlllillissiiiiil" 



e3BSgg3g*g*S**ft** 



i H H i i i i i ii i i i i i i i i i i i i i i i j >>|'|!y ;! 

l|||i|ii||i|fiili MMUHS § 




He 
C 
N 

Ne 
Mg 
Fe 
Ni 



3 4 5 6 7 8 

log (t/yr) 

Figure 4. An evolutionary sequence of a model with neither 
mass loss nor turbulence. Top panel: time evolution of pulsation 
periods. Red dots indicate unstable pulsation modes while green 
crosses are stable modes. The shaded areas are the period ranges 
for the observed p- (lower) and g-modes (upper). Bottom panel: 
time evolution of surface abundances. The top x-axis displays the 
T e ff corresponding to the stellar age on the bottom x-axis. 



is terminated after 10 7 yr because of poor convergence and 
hence long computing time. 

In the upper panels of Fig. [5] we show the effect of these 
mass-loss rates on the stability of the pulsation modes. We 
find that, as expected, for M ^ lO'^Mgyr" 1 pulsations 
cannot be driven. For M = 5 X 10~ 15 Mq yr -1 pulsations are 
driven only during the first 10 5 yr of evolution. This time- 
scale is 100 x smaller than that found by Fontaine et al.| 
( 2006[ ) for a comparable mass- loss rate. It is too short com- 
pared to the evolutionary timescale of 10 s yr to explain the 
fraction of sdBs that are variable (about 10 per cent) as sug- 



gested by Fontaine et al. (20061. The difference is caused by 



their assumption that an equilibrium between gravitational 
settling and radiative levitation can be reached before the 
onset of mass loss. Furthermore, they did not include diffu- 
sion during the evolution while mass is lost. In our models, 
however, diffusion operates at the same time as mass loss and 
both processes start at the ZAHB. We find that, under these 
circumstances, the rates should be M ^ 10~ 15 Moyr _1 if 
pulsations are to be driven for a significant amount of time. 

Remember that, if mass loss were to explain the ob- 
served He abundances, then the r ates should be in the range 
10" 14 < Af/Moyr -1 < 10" 13 ^Fontaine fe Chayer|[l997 



(lower panels) where the surface abundances of He, Fe and 
Ni are plotted. Our results show that these mass-loss rates 
are not consistent with the observed pulsations in sdB stars. 
One could argue that the constant sdBs are experiencing 
mass loss while the pulsators are not (or much less). But 
then there should be a negative correlation between the He 
abundance and the presence of pulsations. This has not been 
observed despite many efforts (e.g. O'Toole fe Heber||2006| 
Blanchette et al.|[2008 1. 

Another possibility is a variable mass loss during stel- 
lar evolution. When the mass-loss rate is low, diffusion can 
build up enough Fe and Ni to drive pulsations. When the 
rate increases, the Fe and Ni reservoir is emptied and the 
star becomes constant. We use the mass loss recipe by |Vink| 
& Cassisi (20021) to estimate how much the mass-loss rate 



could vary. Applying their equation (5) to our models, we 
find that the rate increases gradually from M ~ 4 x 10~ 13 
to 2 x 10" 12 M yr" 1 during the sdB lifetime. Although the 
rates themselves could be overestimated due to an under- 
estimation of the terminal velocity (Unglaub 20081, the de- 
pendence on the stellar parameters (T e ff, L t , M*, Z*) should 
hold ( Vink et al.| [2000). Thus we expect the increase in 
the mass-loss rate, which is mostly due to the increasing 
luminosity as the star evolve^] to remain approximately 
valid. This factor 5 over the sdB lifetime is not sufficient 
to obtain the rates needed to keep He from settling (M > 
10" 14 Mq yr" 1 ) if starting with a rate low enough to drive 
pulsations (M ^ 10~ 15 Mq yr" 1 ). Hence we find it unlikely 
that mass loss is the dominant process for slowing down 
atomic diffusion in sdB stars, although a more sophisticated 
understanding of mass loss is required to draw any definite 
conclusions. 

Fig. [5] also shows that, in the presence of mass loss, the 
surface abundances are poor indicators of whether a star is 
pulsating. For example, during the first 10 6 yr of the simu- 
lation with M — 10~ 14 Mq yr" 1 , both Fe and Ni reach over- 
abundances above 1 dex while no modes are excited. This is 
because, although the surface is enriched, Fe en Ni are actu- 
ally depleted in the driving region^Jln Fig.[6]we plot the time 
evolution of the interior abundances of Fe and Ni. Basically 
what happens is that, as the outermost layers of the star 
are removed, the regions underneath become exposed. Thus 
as time progresses and more mass is peeled away, the sur- 
face abundances are determined by processes that occurred 
deeper in the star. After 10 5 yr the surface is enriched in Fe 
and Ni because the region where accumulation took place 
earlier (AM ~ M At — 10~ 9 Mq) is now exposed. However, 
the driving region at this time corresponds to a region that 
was previously deeper inside the star and from which Fe and 



The dependence of M on Z„ has a small negative effect dur- 
ing the evolution because the metallicity decreases as the lighter 
metals sink. Although one might expect that the mass-loss rate 
increases with Fe/H], one should not forget that the relative role 
of Fe, compared to lighter elements, diminishes for weaker winds 
| |Vink et al.|200l"| l. A detailed study that evaluates the contribu- 
tion of different metals to the mass loss is needed but beyond the 
scope of this work. 

1 The driving is caused by the iron-group opacity bump at 
log(T/K) « 5.3. This corresponds to AM ta lCr lo M in our 
stellar models. 
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Figure 5. Models with mass loss. Upper panels show pulsations as a function of time for mass-loss rates (a) 10 15 , (b) 5 X 10 15 , (c) 
10~ 14 and (d) 10" 13 Moyr" 1 . The red dots are unstable modes while the green crosses are stable. Lower panels show surface abundances 
of He, Fe, and Ni corresponding to these mass-loss rates. The kink in the surface abundances for M = 10 — 15 Mg yr — 1 at 10 4 yr is caused 
by the overshoot of the iron-group convection zone assumed in this particular simulation. 
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Figure 6. Interior abundance profiles of Fe and Ni at different 
stellar ages for the simulation with M = lO _14 M0yr _1 . Lines 
are as in Fig. [3] 



Ni have been transported away. After 10 7 yr so much mass 
is stripped away that the depleted regions are now at the 
surface, causing underabundances of Fe and Ni. The actual 
situation is more complicated because diffusion operates at 
the same time as mass is lost and the diffusion velocities 
change as the star evolves. 

3.3 Atomic diffusion and turbulent mixing 

If mass loss is not responsible for retarding atomic diffusion 



Michaud et al 



then perhaps turbulence is 
that their turbulence model (Eq. with C 
and A Mo = 10" 7 ' 5 



(2011) showed 
= 10 4 , n = 4 
| can explain most observed abun- 
dance anomalies in field and globular cluster sdBs. In their 



t This results in the outer AM 10 — 7 Mq of the star being 
completely mixed. The mixed mass varies slightly with atomic 
species. 



models Fe is near solar throughout the entire mixed re- 
gion, including the driving region. This raises doubt as to 
whether pulsations can be excited. As they conclude them- 
selves, their work has yet to be tested with asteroseismology. 
For this purpose, we run simulations with AMq = 10 -8 ' 5 , 
10 -7 ' 5 and 1O _6 ' 5 M , corresponding to efficient mixing of 
the outer AM w 10" 8 , 10" 7 and 10" 6 M , respectively. The 
other parameters are kept at C = 10 4 and n = 4. 

In Fig. [TJi-c we show the effect of these three turbu- 
lence models on the stability of the pulsation modes. We 
see that mixing the outer mass of AAf w 10 -8 M allows 
some modes to be driven. Although Fe stays near solar, Ni 
can still accumulate and facilitate mode excitation. How- 
ever, the He surface abundance still drops too quickly to 
explain the observations. Mixing more mass (AM « 10 -7 
M ) slows down He settling further but it also prevents the 
accumulation of both Fe and Ni and hence mode excitation. 
Interestingly, if mixing occurs down to AM w 10 -6 M a 
few pulsation modes can be excited after 8 x 10 7 yr. This can 
be understood from Fig. [3] where we see a second, very small 
accumulation of Fe and Ni appearing around AAf ~ 10 -6 
M ) at late stellar ages. However, the number of unstable 
modes is very small and cannot explain the many (dozens) 
frequencies observed in sdB stars. 

It would at first sight seem that turbulence cannot be 
reconciled with mode excitation either. However, this can be 
solved by adjusting the efficiency of turbulence because the 
amount of mixing can vary with atomic species. After all, 
the diffusion velocities of Fe and Ni are much larger than 
for He, so a weak amount of turbulence could retard He set- 
tling while having only a small effect on Fe and Ni levitation. 
This possibility is investigated by varying the parameters C, 
n and AM in Eq. Indeed, we find for C = 100, n = 1 
and AM = 10 -8 ' 2 M that He settles from 0.1 to 10 -4 dur- 
ing the sdB lifetime while Fe and Ni can accumulate enough 
to excite pulsations (see Fig. |7a). Efficient mixing only oc- 
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Figure 7. Models with turbulent mixing. Upper panels show pulsations as a function of time for turbulence models (a) AMq 
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10 4 , n = 4 (b) AMo = 10~ 7 ' 5 M Q , C = 10 4 , n = 4, (c) AM = lO- 6 - 5 M , C = 10 4 , n = 4 and (d) A Mo = 1O- S - 2 M , C = 10 2 , 
n = l. Lower panels show surface abundances of He, Fe, and Ni corresponding to these turbulence models. Symbols and lines are as in 
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Figure 8. Interior abundance profiles of Fe and Ni at different 



stellar ages for the simulation with AMo 
n = 1. Lines are as in Fig. [3] 
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curs in the outer AM w 10 -9 M@, but weak mixing occurs 
up to AM m 10 -6 Mq. This can be seen by the flattening of 
the abundance profiles plotted in Fig. [8] Of course, it can- 
not be guaranteed that the consistency between abundance 
predictions and observations found by |Michaud et aL (2011 1 
can be recovered. However, their turbulence model was fine- 
tuned to keep the Fe abundance near solar, and this remains 
true in our model. 

In Table □ we summarize our results in terms of the 
fraction of its lifetime the sdB star has (i) He surface abun- 



dances in the typically observed range of 10 



0.1, (ii) 



unstable pulsation modes or (iii) both. We see that, of all 
the scenarios we examined, only the sdB models with mix- 
ing in the outer 10 -6 Mq spend a significant time fulfilling 
both conditions. 



4 SUMMARY AND DISCUSSION 

We have constructed sdB evolutionary and seismic models 
with atomic diffusion, including the often neglected radiative 
forces. Our aim is to discriminate between mass loss and tur- 
bulence in sdB stars using non-adiabatic asteroseismology. 
We performed a stability analysis on stellar models with var- 
ious mass-loss rates and turbulence models. We found that 
the mass-loss rates required to match the observed He abun- 
dances are not consistent with observed pulsations. However, 
weak turbulent mixing of the outer 10 -6 Mq can also explain 
the He abundances while still allowing pulsation modes to 
be driven. The presence of unstable modes is very sensitive 
to the turbulence model. This could explain why some sdBs 
are pulsating while others, with similar T e s and logp, are 
constant. 

Although our results favour turbulence over mass loss, 
we should remain cautious because the turbulence is not yet 
supported by a physical model. Possibly, thermohaline mix- 
ing in the presence of a mean molecular weight inversion 



could play a role (Theado et al. 2009). But a detailed in- 



vestigation including the effect of radiative accelerations on 
the /i-inversion instability is still missing. Also, the possibil- 
ity of overshooting beyond the surface convection zones of 
sdB stars should be examined with convection models that 
are more realistic than mixing-length theory. Such studies 
for A stars show that overshoot can cause the regions be- 



tween surface convection zones to completely mix (Kupka 
fe Montgoniery1p002l |Freytag fe Steffen||20u4) . Our mod- 
els based on the mixing-length theory develop a relatively 
broad iron-group convection zone on a very short time-scale. 
Another candidate for causing turbulent mixing could be 
rotation. Although sdB stars typically have slow surface ro- 
tation (usini < lOkrns -1 ), they might hide a rapidly ro- 
tating core, a relic from their previous evolution ( jKawalei] 
fe Hostler|2 005). The differential rotation could cause shear 
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Table 1. Percentage of EHB lifetime (~10 8 yr) spent with He surface abundance in the observed 
range and/or with pulsation modes excited. 



Simulation 


MOTO 


Ml 


M2 


M3 


M4 


Tl 


T2 


T3 


Tl 


(i) He condition 


< 1 


4 


5 


6 


100 


7 


50 


100 


65 


(ii) pulsation condition 


> 99 




< 1 


< 1 





90 





40 


> 99 


(iii) both conditions 





4 


< 1 


< 1 





5 





40 


64 



Simulation (MOTO) is with neither mass loss nor turbulence. The simulations with mass loss are 
(Ml) M = 10- 15 M 9 yr 1 , (M2) M = 5 X 10~ 15 Mq yr" 1 , (M3) M = 10- 14 M e yr- 1 and (M4) 
M = lO- 13 M yr- 1 . The simulations with turbulence are (Tl) AM = 1O~ S ' 5 M , C = 10 4 , 
n = 4, (T2) A Mo = 10" 7 ' 5 M Q , C = 10 4 , n = 4, (T3) A Mo = 10" 6 ' 5 Mq, C = 10 4 , n = 4 and 
(T4) AM = 1O~ 8 - 2 M , C = 10 2 , n = 1. The numbers give the percentage of its life the sdB 
star spends with (i) He surface abundance in the range 10 — 4 — 0.1, (ii) unstable pulsation modes 
or (iii) both conditions fulfilled. Note that simulation (Ml) was terminated after 10 7 yr owing to 
poor convergence, so we only provide an estimate here. 



turbulence but this has never been examined in sdB stars 
as far as we know. We intend to study the influence of these 
additional mixing processes in a forthcoming paper. 

In this study we used observed He abundances to con- 
strain mass-loss rates and turbulent mixing. Metals are not 
trustworthy in this respect because of the possibility of weak 
metallic winds (M < lO _16 M yr _1 ) in sdB stars as sug- 



gested by Unglaub (20081. Such low mass-loss rates could 



still act to change the atmospheric metal abundances with- 
out interfering with mode driving. This complicates or even 
prevents the discrimination between variable and constant 
sdBs on the basis of atmospheric abundances, in agreement 



with hitherto unsuccessful attempts (O 'Toole et al. 2004 



O'Toole fe Heber|2006[ |Blanchette et al.|2008 l 

It is noteworthy that Ni is as important for mode driv- 
ing as Fe. Indeed, the diffusive equilibrium abundances in 
the driving region are comparable despite nickel's lower ini- 
tial abundance. Thus Ni plays an even more important role 
than previous studies of |Jeffery fc Saio (2006J and |Hu et al. 
( 2009 1 indicated because there it was assumed that Ni was 



enhanced by the same factor as Fe. Our models show ex- 
citation of low-degree [t = 1,2) gr-modes at relatively high 
T e ff (29,000 K), so it is likely that the blue-edge problem 
can be resolved. A detailed study of the instability strips is 
underway (Hu et al. in preparation). 
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APPENDIX A: SOLVING THE BURGERS' EQUATIONS 



With the approach of |Thoul et al. (19941, the Burgers' equations without radiative forces, together with the constraints of 
current neutrality and mass conservation, can be expressed as a matrix equation 



k \ 



p / dlnp dlnT v— * dlnC, 
" 'c*i-^ + i>i — h^,- 



dr 



dr 



, 7*j - 



dr 



2S+2 



Most quantities are defined as by Thoul et al. (19941 except for the coefficients Yy and Ay which we define as 

Si 



Hi 



c 



c, 



c 



for i 



Our fij now takes into account the He concentration gradient which was unnecessarily eliminated by Thoul et al. ( 1994 1. The 
coefficients Ay are modified to make use of Paquette et al. ( 1986 1's resistance coefficients (Kij , 2y , 2y and Zy) derived from 
a screened Coulomb potential. For i ■ 



Ay = { 



k^i 
Kij 

^ ^ %ikKi}z'£ik 



for j = i, 

for j = 1,...,S A j ^i, 
for j — i + S, 



For i = 5 + 1, 



A, 



,25, 



Zij-sKi,j-sVi,j-s for .7 = 5+1, 25 A j / i + 5. 

for j = i — 5, 



2.5zi_s ifc Ki_s,fc2:i-s,fc 

— 2.5Zi-S,jKi-S,jXi-S,j 

^2[Ki-s : kyi-s,k(0.8z"_ Syk Xi- S ,k + Yi_ s ,k) - 0-4z"-s,i-sKi-s,i- 
(3 + zl_sj_ s - 0.8z"_ s ,j-s) K i-s,j-syi-s,j-sXi-s,j-s 



ioTj = l,...,SAj^i-S, 
for j = i, 

for j = 5 + l,...,25Aj7H 



with Kij = Kij/Ko, Xij = mj/(mi + rrij), yij = rrn/[rrn + nii) and Yij = 3j/y + z'ijXijViij /mi. Note that the approximations 
( 19941 are retrieved with Zij — 0.6, = 1.3 and z'/j = 0.4 as estimated from a pure Coulomb potential. 



of 



Thoul et al. 



The above modifications were already made by |Hu et al.| ( |2010"| ) but not explicitly described then. In this work we adapted 
Thoul et al. (1994)'s routine to include radiative forces. After some algebra we find that Eqs |2|-([5| can be written as the 
matrix equation 



P 

K 



aimiq la di dlnp dlnT \ -> d In C, 
' +Cti — h Vi — ; — h y j 7,- ■ 



dr 



dr 



3=1 



dr 



2S+2 



(Al) 



The terms on the left-hand side represent the contributions to the diffusion velocity by radiative levitation, gravitational 
settling, thermal diffusion and concentration diffusion. The matrix equation is solved by LU decomposition for each of these 
terms separately and the solutions are combined to give to the total diffusion velocity. In principle, we could also add the 
terms on the left-hand side of Eq. ( Al I beforehand and solve the matrix equation in one step. However, it can be insightful 
to quantify the different contributions to the diffusion velocities. 



